//============================================================================//
//                                                                            //
// TetGen                                                                     //
//                                                                            //
// A Quality Tetrahedral Mesh Generator and A 3D Delaunay Triangulator        //
//                                                                            //
// Version 1.6.0                                                              //
// August 31, 2020                                                            //
//                                                                            //
// Copyright (C) 2002--2020                                                   //
//                                                                            //
// Hang Si                                                                    //
// Research Group: Numerical Mathematics and Scientific Computing             //
// Weierstrass Institute for Applied Analysis and Stochastics (WIAS)          //
// Mohrenstr. 39, 10117 Berlin, Germany                                       //
// si@wias-berlin.de                                                          //
//                                                                            //
// TetGen is a tetrahedral mesh generator. It creates 3d triangulations of    //
//   polyhedral domains. It generates meshes with well-shaped elements whose  //
//   sizes are adapted to the geometric features or user-provided sizing      //
//   functions. It has applications in various applications in scientific     //
//   computing, such as computer graphics (CG), computer-aided design (CAD),  //
//   geometry processing (parametrizations and computer animation), and       //
//   physical simulations (finite element analysis).                          //
//                                                                            //
// TetGen computes (weighted) Delaunay triangulations for three-dimensional   //
//   (weighted) point sets, and constrained Delaunay triangulations for       //
//   three-dimensional polyhedral domains. In the latter case, input edges    //
//   and triangles can be completely preserved in the output meshes.  TetGen  //
//   can refine or coarsen an existing mesh to result in good quality and     //
//   size-adapted mesh according to the geometric features and user-defined   //
//   mesh sizing functions.                                                   //
//                                                                            //
// TetGen implements theoretically proven algorithms for computing the        //
//   Delaunay and constrained Delaunay tetrahedralizations. TetGen achieves   //
//   robustness and efficiency by using advanced techniques in computational  //
//   geometry. A technical paper describes the algorithms and methods         //
//   implemented in TetGen is available in ACM-TOMS, Hang Si ``TetGen, a      //
//   Delaunay-Based Quality Tetrahedral Mesh Generator", ACM Transactions on  //
//   Mathematical Software, February 2015, https://doi.org/10.1145/2629697.   //
//                                                                            //
// TetGen is freely available through the website: http://www.tetgen.org.     //
//   It may be copied, modified, and redistributed for non-commercial use.    //
//   Please consult the file LICENSE for the detailed copyright notices.      //
//============================================================================//
#include <fstream>
#include <iostream>
#include <vector>
#include "tetgenio.h"
using namespace std;

// load tetr. mesh for refinement
bool tetgenio::LoadTetrMesh(string& sName)
{
    // Set initial flags.
    iMeshDim = 3;
    NumCorner = 4;
    NumPointAttribute = 0;  // no point attribute.
    NumTetrAttribute = 1;   // Default no attribute.
    
    iUseIndx = 0;
    string sFile = sName + ".tet";
    ifstream fp(sFile);
    if (!fp.is_open()) {
        cout << "Error : Cannot access file-" << sFile << endl;
        return false;
    }
    cout << "Info : opening the mesh file." << endl;

    int iNumBoud, n;
    fp >> NumPoint >> NumTetr >> iNumBoud;
    if (NumPoint <= 0 || NumTetr <= 0) {
        cout << "Error : Invalid number of points or tetrahedra." << endl;
        fp.close();
        return false;
    }
    pPointList = new REAL[NumPoint * 3];
    // Read the point section.
    REAL x, y, z;
    for (int i = 0; i < NumPoint; ++i) {
        n = 3 * i;
        fp >> x >> y >> z;
        pPointList[n] = x;
        pPointList[n + 1] = y;
        pPointList[n + 2] = z;
    }
    // Read the element
    pTetrList = new int[NumTetr * NumCorner];  // Allocate memory for tetrahedra.
    pTetrAttributeList = new REAL[NumTetr];
    for (int i = 0; i < NumTetr; i++) {
        n = 4 * i;
        fp >> pTetrList[n] >> pTetrList[n + 1] >> pTetrList[n + 2] >> pTetrList[n + 3];
        fp >> pTetrAttributeList[i];
    }
    // Read the boundary
    vector<vector<int>> vTria;
    vector<int> vBoud;
    auto tria = vector<int>(3);
    int id, iNumTria;
    for (int i = 0; i < iNumBoud; ++i) {
        fp >> sFile >> id >> iNumTria;
        for (int j = 0; j < iNumTria; ++j) {
            fp >> tria[0] >> tria[1] >> tria[2];
            vTria.push_back(tria);
            vBoud.push_back(id);
        }
    }
    fp.close();

    if (vTria.empty()) return true;

    NumTriaFace = int(vTria.size());
    trifacelist = new int[NumTriaFace * 3];
    if (trifacelist == (int*)NULL) {
        cout << "Error : Out of memory." << endl;
        exit(1);
    }
    trifacemarkerlist = new int[NumTriaFace];
    if (trifacemarkerlist == (int*)NULL) {
        cout << "Error : Out of memory." << endl;
        exit(1);
    }

    // Read the list of faces.
    for (int i = 0; i < NumTriaFace; ++i) {
        auto p = trifacelist + 3 * i;
        p[0] = vTria[i][0];
        p[1] = vTria[i][1];
        p[2] = vTria[i][2];
        trifacemarkerlist[i] = vBoud[i];
    }
    return true;
}

//============================================================================//
// load_node_call()    Read a list of points from a file.                     //
// 'infile' is the file handle contains the node list.  It may point to a     //
// .node, or .poly or .smesh file. 'markers' indicates each node contains an  //
// additional marker (integer) or not. 'uvflag' indicates each node contains  //
// u,v coordinates or not. It is reuqired by a PSC. 'infilename' is the name  //
// of the file being read,  it is only used in error messages.                //
//                                                                            //
// The 'FirstNumber' (0 or 1) is automatically determined by the number of    //
// the first index of the first point.                                        //
//============================================================================//
bool tetgenio::load_node_call(FILE* infile, int markers, int uvflag, char* infilename)
{
  char inputline[INPUTLINESIZE];
  char *stringptr;
  REAL x, y, z, attrib;
  int firstnode, currentmarker;
  int i, j;

  // Initialize 'pPointList', 'pPointAttributeList', and 'pPointMarkerList'.
  pPointList = new REAL[NumPoint * 3];
  if (pPointList == (REAL *) NULL) {
    printf("Error:  Out of memory.\n");
    exit(1);
  }
  if (NumPointAttribute > 0) {
    pPointAttributeList = new REAL[NumPoint * NumPointAttribute];
    if (pPointAttributeList == (REAL *) NULL) {
        printf("Error:  Out of memory.\n");
        exit(1);
    }
  }
  if (markers) {
    pPointMarkerList = new int[NumPoint];
    if (pPointMarkerList == (int *) NULL) {
        printf("Error:  Out of memory.\n");
        exit(1);
    }
  }
  if (uvflag) {
    pPointParamList = new pointparam[NumPoint];
    if (pPointParamList == NULL) {
        printf("Error:  Out of memory.\n");
        exit(1);
    }
  }

  // Read the point section.
  int index = 0;
  int attribindex = 0;
  for (i = 0; i < NumPoint; i++) {
    stringptr = readnumberline(inputline, infile, infilename);
    if (iUseIndx) {
      if (i == 0) {
        firstnode = (int) strtol (stringptr, &stringptr, 0);
        if ((firstnode == 0) || (firstnode == 1)) {
          FirstNumber = firstnode;
        }
      }
      stringptr = findnextnumber(stringptr);
    } // if (iUseIndx)
    if (*stringptr == '\0') {
      printf("Error:  Point %d has no x coordinate.\n", FirstNumber + i);
      break;
    }
    x = (REAL) strtod(stringptr, &stringptr);
    stringptr = findnextnumber(stringptr);
    if (*stringptr == '\0') {
      printf("Error:  Point %d has no y coordinate.\n", FirstNumber + i);
      break;
    }
    y = (REAL) strtod(stringptr, &stringptr);
    if (iMeshDim == 3) {
      stringptr = findnextnumber(stringptr);
      if (*stringptr == '\0') {
        printf("Error:  Point %d has no z coordinate.\n", FirstNumber + i);
        break;
      }
      z = (REAL) strtod(stringptr, &stringptr);
    } else {
      z = 0.0; // iMeshDim == 2;
    }
    pPointList[index++] = x;
    pPointList[index++] = y;
    pPointList[index++] = z;
    // Read the point attributes.
    for (j = 0; j < NumPointAttribute; j++) {
      stringptr = findnextnumber(stringptr);
      if (*stringptr == '\0') {
        attrib = 0.0;
      } else {
        attrib = (REAL) strtod(stringptr, &stringptr);
      }
      pPointAttributeList[attribindex++] = attrib;
    }
    if (markers) {
      // Read a point marker.
      stringptr = findnextnumber(stringptr);
      if (*stringptr == '\0') {
        currentmarker = 0;
      } else {
        currentmarker = (int) strtol (stringptr, &stringptr, 0);
      }
      pPointMarkerList[i] = currentmarker;
    }
  }
  if (i < NumPoint) {
    // Failed to read points due to some error.
    delete [] pPointList;
    pPointList = (REAL *) NULL;
    if (markers) {
      delete [] pPointMarkerList;
      pPointMarkerList = (int *) NULL;
    }
    if (NumPointAttribute > 0) {
      delete [] pPointAttributeList;
      pPointAttributeList = (REAL *) NULL;
    }
    if (uvflag) {
      delete [] pPointParamList;
      pPointParamList = NULL;
    }
    NumPoint = 0;
    return false;
  }
  return true;
}

//============================================================================//
// load_node()    Load a list of points from a .node file.                    //
//============================================================================//
bool tetgenio::load_node(char* filebasename)
{
  char innodefilename[FILENAMESIZE];
  char inputline[INPUTLINESIZE];

  // Assembling the actual file names we want to open.
  strcpy(innodefilename, filebasename);
  strcat(innodefilename, ".node");

  // Try to open a .node file.
  FILE *infile = fopen(innodefilename, "r");
  if (infile == (FILE *) NULL) {
    printf("  Cannot access file %s.\n", innodefilename);
    return false;
  }
  printf("Opening %s.\n", innodefilename); 

  // Set initial flags.
  iMeshDim = 3;
  NumPointAttribute = 0;  // no point attribute.
  int markers = 0;  // no boundary marker.
  int uvflag = 0; // no uv parameters (required by a PSC). 

  // Read the first line of the file.
  char *stringptr = readnumberline(inputline, infile, innodefilename);
  // Does this file contain an index column?
  stringptr = strstr(inputline, "rbox");
  if (stringptr == NULL) {
      // Read number of points, number of dimensions, number of point
      //   attributes, and number of boundary markers. 
      stringptr = inputline;
      NumPoint = (int)strtol(stringptr, &stringptr, 0);
      stringptr = findnextnumber(stringptr);
      if (*stringptr != '\0') {
        iMeshDim = (int)strtol(stringptr, &stringptr, 0);
      }
      stringptr = findnextnumber(stringptr);
      if (*stringptr != '\0') {
        NumPointAttribute = (int)strtol(stringptr, &stringptr, 0);
      }
      stringptr = findnextnumber(stringptr);
      if (*stringptr != '\0') {
        markers = (int)strtol(stringptr, &stringptr, 0);
      }
      stringptr = findnextnumber(stringptr);
      if (*stringptr != '\0') {
        uvflag = (int)strtol(stringptr, &stringptr, 0);
      }
  } else {
      // It is a rbox (qhull) input file.
      stringptr = inputline;
      // Get the dimension.
      iMeshDim = (int)strtol(stringptr, &stringptr, 0);
      // Get the number of points.
      stringptr = readnumberline(inputline, infile, innodefilename);
      NumPoint = (int)strtol(stringptr, &stringptr, 0);
      // There is no index column.
      iUseIndx = 0;
  }

  // Load the list of nodes.
  bool okflag = load_node_call(infile, markers, uvflag, innodefilename);

  fclose(infile);
  return okflag;
}

//============================================================================//
// load_edge()    Load a list of edges from a .edge file.                     //
//============================================================================//
bool tetgenio::load_edge(char* filebasename)
{
  FILE *infile;
  char inedgefilename[FILENAMESIZE];
  char inputline[INPUTLINESIZE];
  char *stringptr;
  int markers, corner;
  int index;
  int i, j;

  strcpy(inedgefilename, filebasename);
  strcat(inedgefilename, ".edge");

  infile = fopen(inedgefilename, "r");
  if (infile != (FILE *) NULL) {
    printf("Opening %s.\n", inedgefilename);
  } else {
    //printf("  Cannot access file %s.\n", inedgefilename);
    return false;
  }

  // Read number of boundary edges.
  stringptr = readnumberline(inputline, infile, inedgefilename);
  NumEdge = (int) strtol (stringptr, &stringptr, 0);
  if (NumEdge > 0) {
    edgelist = new int[NumEdge * 2];
    if (edgelist == (int *) NULL) {
        printf("Error:  Out of memory.\n");
        exit(1);
    }
    stringptr = findnextnumber(stringptr);
    if (*stringptr == '\0') {
      markers = 0;  // Default value.
    } else {
      markers = (int) strtol (stringptr, &stringptr, 0);
    }
    if (markers > 0) {
      edgemarkerlist = new int[NumEdge];
    }
  }

  // Read the list of edges.
  index = 0;
  for (i = 0; i < NumEdge; i++) {
    // Read edge index and the edge's two endpoints.
    stringptr = readnumberline(inputline, infile, inedgefilename);
    for (j = 0; j < 2; j++) {
      stringptr = findnextnumber(stringptr);
      if (*stringptr == '\0') {
        printf("Error:  Edge %d is missing vertex %d in %s.\n",
               i + FirstNumber, j + 1, inedgefilename);
        printf("Error:  Out of memory.\n");
        exit(1);
      }
      corner = (int) strtol(stringptr, &stringptr, 0);
      if (corner < FirstNumber || corner >= NumPoint + FirstNumber) {
        printf("Error:  Edge %d has an invalid vertex index.\n",
               i + FirstNumber);
        printf("Error:  Out of memory.\n");
        exit(1);
      }
      edgelist[index++] = corner;
    }
    if (NumCorner == 10) {
      // Skip an extra vertex (generated by a previous -o2 option).
      stringptr = findnextnumber(stringptr);
    }
    // Read the edge marker if it has.
    if (markers) {
      stringptr = findnextnumber(stringptr);
      edgemarkerlist[i] = (int) strtol(stringptr, &stringptr, 0);
    }
  }

  fclose(infile);
  return true;
}

//============================================================================//
// load_face()    Load a list of faces (triangles) from a .face file.         //
//============================================================================//
bool tetgenio::load_face(char* filebasename)
{
    string sName = string(filebasename) + ".face";
    // Try to open a .elem file.
    ifstream fp(sName);
    if (!fp.is_open()) {
        return false;
    } else {
        cout << "Info : Opening file-" << sName << endl;
    }

    // Read number of faces, boundary markers.
    int markers;
    fp >> NumTriaFace >> markers;
    if (NumTriaFace > 0) {
        trifacelist = new int[NumTriaFace * 3];
        if (trifacelist == (int *) NULL) {
            cout << "Error : Out of memory." << endl;
            exit(1);
        }
        if (markers > 0) {
            trifacemarkerlist = new int[NumTriaFace];
            if (trifacemarkerlist == (int *) NULL) {
                cout << "Error : Out of memory." << endl;
                exit(1);
            }
        }
    }

    // Read the list of faces.
    int index = 0;
    for (int i = 0; i < NumTriaFace; ++i) {
        auto p = trifacelist + 3 * i;
        fp >> p[0] >> p[1] >> p[2];
        // Read the boundary marker if it exists.
        if (markers > 0) {
            fp >> trifacemarkerlist[i];
        }
    }
    fp.close();
    return true;
}

//============================================================================//
// load_tet()    Load a list of tetrahedra from a .ele file.                  //
//============================================================================//
bool tetgenio::load_tet(char* filebasename)
{
  FILE *infile;
  char infilename[FILENAMESIZE];
  char inputline[INPUTLINESIZE];
  char *stringptr;
  REAL attrib;
  int corner;
  int index, attribindex;
  int i, j;

  strcpy(infilename, filebasename);
  strcat(infilename, ".ele");

  infile = fopen(infilename, "r");
  if (infile != (FILE *) NULL) {
    printf("Opening %s.\n", infilename);
  } else {
    return false;
  }

  // Read number of elements, number of corners (4 or 10), number of
  //   element attributes.
  stringptr = readnumberline(inputline, infile, infilename);
  NumTetr = (int) strtol (stringptr, &stringptr, 0);
  if (NumTetr <= 0) {
    printf("Error:  Invalid number of tetrahedra.\n");
    fclose(infile);
    return false;
  }
  stringptr = findnextnumber(stringptr);
  if (*stringptr == '\0') {
    NumCorner = 4;  // Default read 4 nodes per element.
  } else {
    NumCorner = (int) strtol(stringptr, &stringptr, 0);
  }
  stringptr = findnextnumber(stringptr);
  if (*stringptr == '\0') {
    NumTetrAttribute = 0; // Default no attribute.
  } else {
    NumTetrAttribute = (int) strtol(stringptr, &stringptr, 0);
  }
  if (NumCorner != 4 && NumCorner != 10) {
    printf("Error:  Wrong number of corners %d (should be 4 or 10).\n", NumCorner);
    fclose(infile);
    return false;
  }

  // Allocate memory for tetrahedra.
  pTetrList = new int[NumTetr * NumCorner]; 
  if (pTetrList == (int *) NULL) {
      printf("Error:  Out of memory.\n");
      exit(1);
  }
  // Allocate memory for output tetrahedron attributes if necessary.
  if (NumTetrAttribute > 0) {
    pTetrAttributeList = new REAL[NumTetr *
                                        NumTetrAttribute];
    if (pTetrAttributeList == (REAL *) NULL) {
        printf("Error:  Out of memory.\n");
        exit(1);
    }
  }

  // Read the list of tetrahedra.
  index = 0;
  attribindex = 0;
  for (i = 0; i < NumTetr; i++) {
    // Read tetrahedron index and the tetrahedron's corners.
    stringptr = readnumberline(inputline, infile, infilename);
    for (j = 0; j < NumCorner; j++) {
      stringptr = findnextnumber(stringptr);
      if (*stringptr == '\0') {
        printf("Error:  Tetrahedron %d is missing vertex %d in %s.\n",
               i + FirstNumber, j + 1, infilename);
        printf("Error:  Out of memory.\n");
        exit(1);
      }
      corner = (int) strtol(stringptr, &stringptr, 0);
      if (corner < FirstNumber || corner >= NumPoint + FirstNumber) {
        printf("Error:  Tetrahedron %d has an invalid vertex index.\n",
               i + FirstNumber);
        printf("Error:  Out of memory.\n");
        exit(1);
      }
      pTetrList[index++] = corner;
    }
    // Read the tetrahedron's attributes.
    for (j = 0; j < NumTetrAttribute; j++) {
      stringptr = findnextnumber(stringptr);
      if (*stringptr == '\0') {
        attrib = 0.0;
      } else {
        attrib = (REAL) strtod(stringptr, &stringptr);
      }
      pTetrAttributeList[attribindex++] = attrib;
    }
  }

  fclose(infile);

  return true;
}

//============================================================================//
// load_vol()    Load a list of volume constraints from a .vol file.          //
//============================================================================//
bool tetgenio::load_vol(char* filebasename)
{
    string sName = string(filebasename) + ".vol";
    // Try to open a .elem file.
    ifstream fp(sName);
    if (!fp.is_open()) {
        return false;
    } else {
        cout << "Info : Opening file-" << sName << endl;
    }

    // Read number of tetrahedra.
    pTetrVolumeList = new REAL[NumTetr];
    if (pTetrVolumeList == (REAL*)NULL) {
        cout << "Error : Out of memory." << endl;
        exit(1);
    }
    // Read the list of volume constraints.
    for (int i = 0; i < NumTetr; i++) {
        fp >> pTetrVolumeList[i];
    }
    fp.close();
    return true;
}

//============================================================================//
// load_var()    Load constraints applied on facets, segments, and nodes      //
//               from a .var file.                                            //
//============================================================================//
bool tetgenio::load_var(char* filebasename)
{
  FILE *infile;
  char varfilename[FILENAMESIZE];
  char inputline[INPUTLINESIZE];
  char *stringptr;
  int index;
  int i;

  // Variant constraints are saved in file "filename.var".
  strcpy(varfilename, filebasename);
  strcat(varfilename, ".var");
  infile = fopen(varfilename, "r");
  if (infile != (FILE *) NULL) {
    printf("Opening %s.\n", varfilename);
  } else {
    return false;
  }

  // Read the facet constraint section.
  stringptr = readnumberline(inputline, infile, varfilename);
  if (stringptr == NULL) {
    // No region list, return.
    fclose(infile);
    return true;
  }
  if (*stringptr != '\0') {
    NumFacetConstraint = (int) strtol (stringptr, &stringptr, 0);
  } else {
    NumFacetConstraint = 0;
  }
  if (NumFacetConstraint > 0) {
    // Initialize 'pFacetConstraintList'.
    pFacetConstraintList = new REAL[NumFacetConstraint * 2];
    index = 0;
    for (i = 0; i < NumFacetConstraint; i++) {
      stringptr = readnumberline(inputline, infile, varfilename);
      stringptr = findnextnumber(stringptr);
      if (*stringptr == '\0') {
        printf("Error:  facet constraint %d has no facet marker.\n",FirstNumber + i);
        break;
      } else {
        pFacetConstraintList[index++] = (REAL) strtod(stringptr, &stringptr);
      }
      stringptr = findnextnumber(stringptr);
      if (*stringptr == '\0') {
        printf("Error:  facet constraint %d has no maximum area bound.\n",FirstNumber + i);
        break;
      } else {
        pFacetConstraintList[index++] = (REAL) strtod(stringptr, &stringptr);
      }
    }
    if (i < NumFacetConstraint) {
      // This must be caused by an error.
      fclose(infile);
      return false;
    }
  }

  // Read the segment constraint section.
  stringptr = readnumberline(inputline, infile, varfilename);
  if (stringptr == NULL) {
    // No segment list, return.
    fclose(infile);
    return true;
  }
  if (*stringptr != '\0') {
    NumSegmentConstraint = (int) strtol (stringptr, &stringptr, 0);
  } else {
    NumSegmentConstraint = 0;
  }
  if (NumSegmentConstraint > 0) {
    // Initialize 'pSegmentConstraintList'.
    pSegmentConstraintList = new REAL[NumSegmentConstraint * 3];
    index = 0;
    for (i = 0; i < NumSegmentConstraint; i++) {
      stringptr = readnumberline(inputline, infile, varfilename);
      stringptr = findnextnumber(stringptr);
      if (*stringptr == '\0') {
        printf("Error:  segment constraint %d has no frist endpoint.\n",FirstNumber + i);
        break;
      } else {
        pSegmentConstraintList[index++] = (REAL) strtod(stringptr, &stringptr);
      }
      stringptr = findnextnumber(stringptr);
      if (*stringptr == '\0') {
        printf("Error:  segment constraint %d has no second endpoint.\n",FirstNumber + i);
        break;
      } else {
        pSegmentConstraintList[index++] = (REAL) strtod(stringptr, &stringptr);
      }
      stringptr = findnextnumber(stringptr);
      if (*stringptr == '\0') {
        printf("Error:  segment constraint %d has no maximum length bound.\n",FirstNumber + i);
        break;
      } else {
        pSegmentConstraintList[index++] = (REAL) strtod(stringptr, &stringptr);
      }
    }
    if (i < NumSegmentConstraint) {
      // This must be caused by an error.
      fclose(infile);
      return false;
    }
  }

  fclose(infile);
  return true;
}

//============================================================================//
// load_mtr()    Load a size specification map from a .mtr file.              //
//============================================================================//
bool tetgenio::load_mtr(char* filebasename)
{
  FILE *infile;
  char mtrfilename[FILENAMESIZE];
  char inputline[INPUTLINESIZE];
  char *stringptr;
  REAL mtr;
  int ptnum;
  int mtrindex;
  int i, j;

  strcpy(mtrfilename, filebasename);
  strcat(mtrfilename, ".mtr");
  infile = fopen(mtrfilename, "r");
  if (infile != (FILE *) NULL) {
    printf("Opening %s.\n", mtrfilename);
  } else {
    return false;
  }

  // Read the number of points.
  stringptr = readnumberline(inputline, infile, mtrfilename);
  ptnum = (int) strtol (stringptr, &stringptr, 0);
  if (ptnum != NumPoint) {
    printf("  !! Point numbers are not equal. Ignored.\n");
    fclose(infile);
    return false;
  }
  // Read the number of columns (1, 3, or 6).
  stringptr = findnextnumber(stringptr); // Skip number of points.
  if (*stringptr != '\0') {
    numberofpointmtrs = (int) strtol (stringptr, &stringptr, 0);
  }
  if ((numberofpointmtrs != 1) && (numberofpointmtrs != 3) &&
      (numberofpointmtrs != 6)) {
    // Column number doesn't match. 
    numberofpointmtrs = 0;
    printf("  !! Metric size does not match (1, 3, or 6). Ignored.\n");
    fclose(infile);
    return false;
  }

  // Allocate space for pPointMtrList.
  pPointMtrList = new REAL[NumPoint * numberofpointmtrs];
  if (pPointMtrList == (REAL *) NULL) {
      printf("Error:  Out of memory.\n");
      exit(1);
  }
  mtrindex = 0;
  for (i = 0; i < NumPoint; i++) {
      // Read metrics.
      stringptr = readnumberline(inputline, infile, mtrfilename);
      for (j = 0; j < numberofpointmtrs; j++) {
          if (*stringptr == '\0') {
              printf("Error:  Metric %d is missing value #%d in %s.\n",i + FirstNumber, j + 1, mtrfilename);
              printf("Error:  Out of memory.\n");
              exit(1);
          }
          mtr = (REAL) strtod(stringptr, &stringptr);
          pPointMtrList[mtrindex++] = mtr;
          stringptr = findnextnumber(stringptr);
      }
  }
  fclose(infile);
  return true;
}

//============================================================================//
// load_elem()    Load a list of refine elements from an .elem file.          //
//============================================================================//
bool tetgenio::load_elem(char* filebasename)
{
    string sName = string(filebasename) + ".elem";
    // Try to open a .elem file.
    ifstream fp(sName);
    if (!fp.is_open()) {
        return false;
    } else {
        cout << "Info : Opening file-" << sName << endl;
    }

    int elenum = 0;
    float growth_ratio = 0.0;
    fp >> elenum >> growth_ratio;

    if (0 == elenum) {
        fp.close();
        return false;
    }
    refine_elem_list = new int[elenum * 4];
    NumRefineElem = elenum;

    for (int i = 0; i < elenum; ++i) {
      auto idx = refine_elem_list + i * 4;
      fp >> idx[0] >> idx[1] >> idx[2] >> idx[3];
    }
    fp.close();
    return true;
}

//============================================================================//
// load_poly()    Load a PL complex from a .poly or a .smesh file.            //
//============================================================================//
bool tetgenio::load_poly(char* filebasename)
{
  FILE *infile;
  char inpolyfilename[FILENAMESIZE];
  char insmeshfilename[FILENAMESIZE];
  char inputline[INPUTLINESIZE];
  char *stringptr, *infilename;
  int smesh, markers, uvflag, currentmarker;
  int index;
  int i, j, k;

  // Assembling the actual file names we want to open.
  strcpy(inpolyfilename, filebasename);
  strcpy(insmeshfilename, filebasename);
  strcat(inpolyfilename, ".poly");
  strcat(insmeshfilename, ".smesh");

  // First assume it is a .poly file.
  smesh = 0;
  // Try to open a .poly file.
  infile = fopen(inpolyfilename, "r");
  if (infile == (FILE *) NULL) {
    // .poly doesn't exist! Try to open a .smesh file.
    infile = fopen(insmeshfilename, "r");
    if (infile == (FILE *) NULL) {
      printf("  Cannot access file %s and %s.\n",
             inpolyfilename, insmeshfilename);
      return false;
    } else {
      printf("Opening %s.\n", insmeshfilename);
      infilename = insmeshfilename;
    }
    smesh = 1;
  } else {
    printf("Opening %s.\n", inpolyfilename);
    infilename = inpolyfilename;
  }

  // Initialize the default values.
  iMeshDim = 3;  // Three-dimensional coordinates.
  NumPointAttribute = 0;  // no point attribute.
  markers = 0;  // no boundary marker.
  uvflag = 0; // no uv parameters (required by a PSC).

  // Read number of points, number of dimensions, number of point
  //   attributes, and number of boundary markers.
  stringptr = readnumberline(inputline, infile, infilename);
  NumPoint = (int) strtol (stringptr, &stringptr, 0);
  stringptr = findnextnumber(stringptr);
  if (*stringptr != '\0') {
    iMeshDim = (int) strtol (stringptr, &stringptr, 0);      
  }
  stringptr = findnextnumber(stringptr);
  if (*stringptr != '\0') {
    NumPointAttribute = (int) strtol (stringptr, &stringptr, 0);
  }
  stringptr = findnextnumber(stringptr);
  if (*stringptr != '\0') {
    markers = (int) strtol (stringptr, &stringptr, 0);
  }
  if (*stringptr != '\0') {
    uvflag = (int) strtol (stringptr, &stringptr, 0);
  }

  if (NumPoint > 0) {
    // Load the list of nodes.
    if (!load_node_call(infile, markers, uvflag, infilename)) {
      fclose(infile);
      return false;
    }
  } else {
    // If the .poly or .smesh file claims there are zero points, that
    //   means the points should be read from a separate .node file.
    if (!load_node(filebasename)) {
      fclose(infile);
      return false;
    }
  }

  if ((iMeshDim != 3) && (iMeshDim != 2)) {
    printf("Input error:  TetGen only works for 2D & 3D point sets.\n");
    fclose(infile);
    return false;
  }
  if (NumPoint < (iMeshDim + 1)) {
    printf("Input error:  TetGen needs at least %d points.\n", iMeshDim + 1);
    fclose(infile);
    return false;
  }

  facet *f;
  polygon *p;

  if (iMeshDim == 3) {

    // Read number of facets and number of boundary markers.
    stringptr = readnumberline(inputline, infile, infilename);
    if (stringptr == NULL) {
      // No facet list, return.
      fclose(infile);
      return true;
    }
    NumFacet = (int) strtol (stringptr, &stringptr, 0);
    if (NumFacet <= 0) {
      // No facet list, return.
      fclose(infile);
      return true;
    }
    stringptr = findnextnumber(stringptr);
    if (*stringptr == '\0') {
      markers = 0;  // no boundary marker.
    } else {
      markers = (int) strtol (stringptr, &stringptr, 0);
    }

    // Initialize the 'pFacetList', 'pFacetMarkerList'.
    pFacetList = new facet[NumFacet];
    if (markers == 1) {
      pFacetMarkerList = new int[NumFacet];
    }

    // Read data into 'pFacetList', 'pFacetMarkerList'.
    if (smesh == 0) {
      // Facets are in .poly file format.
      for (i = 1; i <= NumFacet; i++) {
        f = &(pFacetList[i - 1]);
        init(f);
        f->NumHole = 0;
        currentmarker = 0;
        // Read number of polygons, number of holes, and a boundary marker.
        stringptr = readnumberline(inputline, infile, infilename);
        f->NumPolygon = (int) strtol (stringptr, &stringptr, 0);
        stringptr = findnextnumber(stringptr);
        if (*stringptr != '\0') {
          f->NumHole = (int) strtol (stringptr, &stringptr, 0);
          if (markers == 1) {
            stringptr = findnextnumber(stringptr);
            if (*stringptr != '\0') {
              currentmarker = (int) strtol(stringptr, &stringptr, 0);
            }
          }
        }
        // Initialize facetmarker if it needs.
        if (markers == 1) {
          pFacetMarkerList[i - 1] = currentmarker; 
        }
        // Each facet should has at least one polygon.
        if (f->NumPolygon <= 0) {
          printf("Error:  Wrong number of polygon in %d facet.\n", i);
          break; 
        }
        // Initialize the 'f->pPolygonList'.
        f->pPolygonList = new polygon[f->NumPolygon];
        // Go through all polygons, read in their vertices.
        for (j = 1; j <= f->NumPolygon; j++) {
          p = &(f->pPolygonList[j - 1]);
          init(p);
          // Read number of vertices of this polygon.
          stringptr = readnumberline(inputline, infile, infilename);
          p->NumVertex = (int) strtol(stringptr, &stringptr, 0);
          if (p->NumVertex < 1) {
            printf("Error:  Wrong polygon %d in facet %d\n", j, i);
            break;
          }
          // Initialize 'p->pVertexList'.
          p->pVertexList = new int[p->NumVertex];
          // Read all vertices of this polygon.
          for (k = 1; k <= p->NumVertex; k++) {
            stringptr = findnextnumber(stringptr);
            if (*stringptr == '\0') {
              // Try to load another non-empty line and continue to read the
              //   rest of vertices.
              stringptr = readnumberline(inputline, infile, infilename);
              if (*stringptr == '\0') {
                printf("Error: Missing %d endpoints of polygon %d in facet %d",
                       p->NumVertex - k, j, i);
                break;
              }
            }
            p->pVertexList[k - 1] = (int) strtol (stringptr, &stringptr, 0);
          }
        } 
        if (j <= f->NumPolygon) {
          // This must be caused by an error. However, there're j - 1
          //   polygons have been read. Reset the 'f->numberofpolygon'.
          if (j == 1) {
            // This is the first polygon.
            delete [] f->pPolygonList;
          }
          f->NumPolygon = j - 1;
          // No hole will be read even it exists.
          f->NumHole = 0;
          break;
        }
        // If this facet has hole pints defined, read them.
        if (f->NumHole > 0) {
          // Initialize 'f->pHoleList'.
          f->pHoleList = new REAL[f->NumHole * 3];
          // Read the holes' coordinates.
          index = 0;
          for (j = 1; j <= f->NumHole; j++) {
            stringptr = readnumberline(inputline, infile, infilename);
            for (k = 1; k <= 3; k++) {
              stringptr = findnextnumber(stringptr);
              if (*stringptr == '\0') {
                printf("Error:  Hole %d in facet %d has no coordinates", j, i);
                break;
              }
              f->pHoleList[index++] = (REAL) strtod (stringptr, &stringptr);
            }
            if (k <= 3) {
              // This must be caused by an error.
              break;
            }
          }
          if (j <= f->NumHole) {
            // This must be caused by an error.
            break;
          }
        }
      }
      if (i <= NumFacet) {
        // This must be caused by an error.
        NumFacet = i - 1;
        fclose(infile);
        return false;
      }
    } else { // poly == 0
      // Read the facets from a .smesh file.
      for (i = 1; i <= NumFacet; i++) {
        f = &(pFacetList[i - 1]);
        init(f);
        // Initialize 'f->pFacetList'. In a .smesh file, each pFacetList only
        //   contains exactly one polygon, no hole.
        f->NumPolygon = 1;
        f->pPolygonList = new polygon[f->NumPolygon];
        p = &(f->pPolygonList[0]);
        init(p);
        // Read number of vertices of this polygon.
        stringptr = readnumberline(inputline, infile, insmeshfilename);
        p->NumVertex = (int) strtol (stringptr, &stringptr, 0);
        if (p->NumVertex < 1) {
          printf("Error:  Wrong number of vertex in facet %d\n", i);
          break;
        }
        // Initialize 'p->pVertexList'.
        p->pVertexList = new int[p->NumVertex];
        for (k = 1; k <= p->NumVertex; k++) {
          stringptr = findnextnumber(stringptr);
          if (*stringptr == '\0') {
            // Try to load another non-empty line and continue to read the
            //   rest of vertices.
            stringptr = readnumberline(inputline, infile, infilename);
            if (*stringptr == '\0') {
              printf("Error:  Missing %d endpoints in facet %d",
                     p->NumVertex - k, i);
              break;
            }
          }
          p->pVertexList[k - 1] = (int) strtol (stringptr, &stringptr, 0);
        }
        if (k <= p->NumVertex) {
          // This must be caused by an error.
          break;
        }
        // Read facet's boundary marker at last.
        if (markers == 1) {
          stringptr = findnextnumber(stringptr);
          if (*stringptr == '\0') {
            currentmarker = 0;
          } else {
            currentmarker = (int) strtol(stringptr, &stringptr, 0);
          }
          pFacetMarkerList[i - 1] = currentmarker;
        }
      }
      if (i <= NumFacet) {
        // This must be caused by an error.
        NumFacet = i - 1;
        fclose(infile);
        return false;
      }
    }

    // Read the hole section.
    stringptr = readnumberline(inputline, infile, infilename);
    if (stringptr == NULL) {
      // No hole list, return.
      fclose(infile);
      return true;
    }
    if (*stringptr != '\0') {
      NumHole = (int) strtol (stringptr, &stringptr, 0);
    } else {
      NumHole = 0;
    }
    if (NumHole > 0) {
      // Initialize 'pHoleList'.
      pHoleList = new REAL[NumHole * 3];
      for (i = 0; i < 3 * NumHole; i += 3) {
        stringptr = readnumberline(inputline, infile, infilename);
        stringptr = findnextnumber(stringptr);
        if (*stringptr == '\0') {
          printf("Error:  Hole %d has no x coord.\n", FirstNumber + (i / 3));
          break;
        } else {
          pHoleList[i] = (REAL) strtod(stringptr, &stringptr);
        }
        stringptr = findnextnumber(stringptr);
        if (*stringptr == '\0') {
          printf("Error:  Hole %d has no y coord.\n", FirstNumber + (i / 3));
          break;
        } else {
          pHoleList[i + 1] = (REAL) strtod(stringptr, &stringptr);
        }
        stringptr = findnextnumber(stringptr);
        if (*stringptr == '\0') {
          printf("Error:  Hole %d has no z coord.\n", FirstNumber + (i / 3));
          break;
        } else {
          pHoleList[i + 2] = (REAL) strtod(stringptr, &stringptr);
        }
      }
      if (i < 3 * NumHole) {
        // This must be caused by an error.
        fclose(infile);
        return false;
      }
    }

    // Read the region section.  The 'region' section is optional, if we
    //   don't reach the end-of-file, try read it in.
    stringptr = readnumberline(inputline, infile, NULL);
    if (stringptr != (char *) NULL && *stringptr != '\0') {
      NumRegion = (int) strtol (stringptr, &stringptr, 0);
    } else {
      NumRegion = 0;
    }
    if (NumRegion > 0) {
      // Initialize 'pRegionList'.
      pRegionList = new REAL[NumRegion * 5];
      index = 0;
      for (i = 0; i < NumRegion; i++) {
        stringptr = readnumberline(inputline, infile, infilename);
        stringptr = findnextnumber(stringptr);
        if (*stringptr == '\0') {
          printf("Error:  Region %d has no x coordinate.\n", FirstNumber + i);
          break;
        } else {
          pRegionList[index++] = (REAL) strtod(stringptr, &stringptr);
        }
        stringptr = findnextnumber(stringptr);
        if (*stringptr == '\0') {
          printf("Error:  Region %d has no y coordinate.\n", FirstNumber + i);
          break;
        } else {
          pRegionList[index++] = (REAL) strtod(stringptr, &stringptr);
        }
        stringptr = findnextnumber(stringptr);
        if (*stringptr == '\0') {
          printf("Error:  Region %d has no z coordinate.\n", FirstNumber + i);
          break;
        } else {
          pRegionList[index++] = (REAL) strtod(stringptr, &stringptr);
        }
        stringptr = findnextnumber(stringptr);
        if (*stringptr == '\0') {
          printf("Error:  Region %d has no region attrib.\n", FirstNumber + i);
          break;
        } else {
          pRegionList[index++] = (REAL) strtod(stringptr, &stringptr);
        }
        stringptr = findnextnumber(stringptr);
        if (*stringptr == '\0') {
          pRegionList[index] = pRegionList[index - 1];
        } else {
          pRegionList[index] = (REAL) strtod(stringptr, &stringptr);
        }
        index++;
      }
      if (i < NumRegion) {
        // This must be caused by an error.
        fclose(infile);
        return false;
      }
    }

  }

  // End of reading poly/smesh file.
  fclose(infile);
  return true;
}

//============================================================================//
//                                                                            //
// load_ply()    Load a polyhedron from a .ply file.                          //
//                                                                            //
// This is a simplified version of reading .ply files, which only reads the   //
// set of vertices and the set of faces. Other informations (such as color,   //
// material, texture, etc) in .ply file are ignored. Complete routines for    //
// reading and writing ,ply files are available from: http://www.cc.gatech.   //
// edu/projects/large_models/ply.html.  Except the header section, ply file   //
// format has exactly the same format for listing vertices and polygons as    //
// off file format.                                                           //
//                                                                            //
//============================================================================//

bool tetgenio::load_ply(char* filebasename)
{
  FILE *fp;
  tetgenio::facet *f;
  tetgenio::polygon *p;
  char infilename[FILENAMESIZE];
  char buffer[INPUTLINESIZE];
  char *bufferp, *str;
  double *coord;
  int endheader = 0, format = 0;
  int nverts = 0, iverts = 0;
  int nfaces = 0, ifaces = 0;
  int line_count = 0, i;

  // Default, the ply file's index is from '0'. We check it by remembering the
  //   smallest index we found in the file. It should be either 0 or 1.
  int smallestidx = 0; 

  strncpy(infilename, filebasename, FILENAMESIZE - 1);
  infilename[FILENAMESIZE - 1] = '\0';
  if (infilename[0] == '\0') {
    printf("Error:  No filename.\n");
    return false;
  }
  if (strcmp(&infilename[strlen(infilename) - 4], ".ply") != 0) {
    strcat(infilename, ".ply");
  }

  if (!(fp = fopen(infilename, "r"))) {
    printf("Error:  Unable to open file %s\n", infilename);
    return false;
  }
  printf("Opening %s.\n", infilename);

  while ((bufferp = readline(buffer, fp, &line_count)) != NULL) {
    if (!endheader) {
      // Find if it is the keyword "end_header".
      str = strstr(bufferp, "end_header");
      // strstr() is case sensitive.
      if (!str) str = strstr(bufferp, "End_header");
      if (!str) str = strstr(bufferp, "End_Header");
      if (str) {
        // This is the end of the header section.
        endheader = 1; 
        continue;
      }
      // Parse the number of vertices and the number of faces.
      if (nverts == 0 || nfaces == 0) {
        // Find if it si the keyword "element".
        str = strstr(bufferp, "element");
        if (!str) str = strstr(bufferp, "Element");
        if (str) {
          bufferp = findnextfield(str);
          if (*bufferp == '\0') {
            printf("Syntax error reading element type on line%d in file %s\n",
                   line_count, infilename);
            fclose(fp);
            return false;
          }
          if (nverts == 0) {
            // Find if it is the keyword "vertex".
            str = strstr(bufferp, "vertex");
            if (!str) str = strstr(bufferp, "Vertex");
            if (str) {
              bufferp = findnextnumber(str);
              if (*bufferp == '\0') {
                printf("Syntax error reading vertex number on line");
                printf(" %d in file %s\n", line_count, infilename);
                fclose(fp);
                return false;
              }
              nverts = (int) strtol(bufferp, &bufferp, 0);
              // Allocate memory for 'tetgenio'
              if (nverts > 0) {
                NumPoint = nverts;
                pPointList = new REAL[nverts * 3];
                smallestidx = nverts + 1; // A big enough index.
              }
            }
          }
          if (nfaces == 0) {
            // Find if it is the keyword "face".
            str = strstr(bufferp, "face");
            if (!str) str = strstr(bufferp, "Face");
            if (str) {
              bufferp = findnextnumber(str);
              if (*bufferp == '\0') {
                printf("Syntax error reading face number on line");
                printf(" %d in file %s\n", line_count, infilename);
                fclose(fp);
                return false;
              }
              nfaces = (int) strtol(bufferp, &bufferp, 0);
              // Allocate memory for 'tetgenio'
              if (nfaces > 0) {        
                NumFacet = nfaces;
                pFacetList = new tetgenio::facet[nfaces];
              }
            }
          }
        } // It is not the string "element". 
      }
      if (format == 0) {
        // Find the keyword "format".
        str = strstr(bufferp, "format");
        if (!str) str = strstr(bufferp, "Format");
        if (str) {
          format = 1;
          bufferp = findnextfield(str);
          // Find if it is the string "ascii".
          str = strstr(bufferp, "ascii");
          if (!str) str = strstr(bufferp, "ASCII");
          if (!str) {
            printf("This routine only reads ascii format of ply files.\n");
            printf("Hint: You can convert the binary to ascii format by\n");
            printf("  using the provided ply tools:\n");
            printf("  ply2ascii < %s > ascii_%s\n", infilename, infilename);
            fclose(fp);
            return false;
          }
        }
      }
    } else if (iverts < nverts) {
      // Read vertex coordinates
      coord = &pPointList[iverts * 3];
      for (i = 0; i < 3; i++) {
        if (*bufferp == '\0') {
          printf("Syntax error reading vertex coords on line %d in file %s\n",
                 line_count, infilename);
          fclose(fp);
          return false;
        }
        coord[i] = (REAL) strtod(bufferp, &bufferp);
        bufferp = findnextnumber(bufferp);
      }
      iverts++;
    } else if (ifaces < nfaces) {
      // Get next face
      f = &pFacetList[ifaces];
      init(f);      
      // In .off format, each facet has one polygon, no hole.
      f->NumPolygon = 1;
      f->pPolygonList = new tetgenio::polygon[1];
      p = &f->pPolygonList[0];
      init(p);
      // Read the number of vertices, it should be greater than 0.
      p->NumVertex = (int) strtol(bufferp, &bufferp, 0);
      if (p->NumVertex == 0) {
        printf("Syntax error reading polygon on line %d in file %s\n", line_count, infilename);
        fclose(fp);
        return false;
      }
      // Allocate memory for face vertices
      p->pVertexList = new int[p->NumVertex];
      for (i = 0; i < p->NumVertex; i++) {
        bufferp = findnextnumber(bufferp);
        if (*bufferp == '\0') {
          printf("Syntax error reading polygon on line %d in file %s\n", line_count, infilename);
          fclose(fp);
          return false;
        }
        p->pVertexList[i] = (int) strtol(bufferp, &bufferp, 0);
        if (p->pVertexList[i] < smallestidx) {
          smallestidx = p->pVertexList[i];
        }
      }
      ifaces++;
    } else {
      // Should never get here
      printf("Found extra text starting at line %d in file %s\n", line_count, infilename);
      break;
    }
  }
  fclose(fp);

  // Decide the FirstNumber of the index.
  if (smallestidx == 0) {
    FirstNumber = 0;  
  } else if (smallestidx == 1) {
    FirstNumber = 1;
  } else {
    printf("A wrong smallest index (%d) was detected in file %s\n", smallestidx, infilename);
    return false;
  }

  if (iverts != nverts) {
    printf("Expected %d vertices, but read only %d vertices in file %s\n", nverts, iverts, infilename);
    return false;
  }
  if (ifaces != nfaces) {
    printf("Expected %d faces, but read only %d faces in file %s\n", nfaces, ifaces, infilename);
    return false;
  }

  return true;
}

//============================================================================//
// load_plc()    Load a piecewise linear complex from file(s).                //
//============================================================================//
bool tetgenio::load_plc(char* filebasename, int object)
{
    bool success;

    if (object == (int) tetgenbehavior::NODES) {
        success = load_node(filebasename);
    } else if (object == (int) tetgenbehavior::POLY) {
        success = load_poly(filebasename);
    } else if (object == (int) tetgenbehavior::PLY) {
        success = load_ply(filebasename);
    } else {
        success = load_poly(filebasename);
    }

    if (success) {
        // Try to load the following files (.edge, .var, .mtr).
        load_edge(filebasename);
        load_var(filebasename);
        load_mtr(filebasename);
    }
    return success;
}

//============================================================================//
// load_mesh()    Load a tetrahedral mesh from file(s).                       //
//============================================================================//
bool tetgenio::load_tetmesh(char* filebasename, int object)
{
    bool success = load_node(filebasename);
    if (success) {
        success = load_tet(filebasename);
    }
    if (success) {
        // Try to load the following files (.face, .edge, .vol).
        load_face(filebasename);
        load_edge(filebasename);
        load_vol(filebasename);

        // Try to load the following files (.var, .mtr).
        load_var(filebasename);
        load_mtr(filebasename);
        load_elem(filebasename);
    }
    return success;
}

//============================================================================//
// save_nodes()    Save points to a .node file.                               //
//============================================================================//
void tetgenio::save_nodes(const char *filebasename)
{
  FILE *fout;
  char outnodefilename[FILENAMESIZE];
  char outmtrfilename[FILENAMESIZE];
  int i, j;

  sprintf(outnodefilename, "%s.node", filebasename);
  printf("Saving nodes to %s\n", outnodefilename);
  fout = fopen(outnodefilename, "w");
  fprintf(fout, "%d  %d  %d  %d\n", NumPoint, iMeshDim,
          NumPointAttribute, pPointMarkerList != NULL ? 1 : 0);
  for (i = 0; i < NumPoint; i++) {
    if (iMeshDim == 2) {
      fprintf(fout, "%d  %.16g  %.16g", i + FirstNumber, pPointList[i * 3],
              pPointList[i * 3 + 1]);
    } else {
      fprintf(fout, "%d  %.16g  %.16g  %.16g", i + FirstNumber,
              pPointList[i * 3], pPointList[i * 3 + 1], pPointList[i * 3 + 2]);
    }
    for (j = 0; j < NumPointAttribute; j++) {
      fprintf(fout, "  %.16g", 
              pPointAttributeList[i * NumPointAttribute + j]);
    }
    if (pPointMarkerList != NULL) {
      fprintf(fout, "  %d", pPointMarkerList[i]);
    }
    fprintf(fout, "\n");
  }
  fclose(fout);

  // If the point metrics exist, output them to a .mtr file.
  if ((numberofpointmtrs > 0) && (pPointMtrList != (REAL *) NULL)) {
    sprintf(outmtrfilename, "%s.mtr", filebasename);
    printf("Saving metrics to %s\n", outmtrfilename);
    fout = fopen(outmtrfilename, "w");
    fprintf(fout, "%d  %d\n", NumPoint, numberofpointmtrs);
    for (i = 0; i < NumPoint; i++) {
      for (j = 0; j < numberofpointmtrs; j++) {
        fprintf(fout, "%.16g ", pPointMtrList[i * numberofpointmtrs + j]);
      }
      fprintf(fout, "\n");
    }
    fclose(fout);
  }
}

//============================================================================//
//                                                                            //
// save_elements()    Save elements to a .ele file.                           //
//                                                                            //
//============================================================================//

void tetgenio::save_elements(const char* filebasename)
{
  FILE *fout;
  char outelefilename[FILENAMESIZE];
  int i, j;

  sprintf(outelefilename, "%s.ele", filebasename);
  printf("Saving elements to %s\n", outelefilename);
  fout = fopen(outelefilename, "w");
  if (iMeshDim == 3) {
    fprintf(fout, "%d  %d  %d\n", NumTetr, NumCorner,
            NumTetrAttribute);
    for (i = 0; i < NumTetr; i++) {
      fprintf(fout, "%d", i + FirstNumber);
      for (j = 0; j < NumCorner; j++) {
        fprintf(fout, "  %5d", pTetrList[i * NumCorner + j]);
      }
      for (j = 0; j < NumTetrAttribute; j++) {
        fprintf(fout, "  %g",
          pTetrAttributeList[i * NumTetrAttribute + j]);
      }
      fprintf(fout, "\n");
    }
  } else {
    // Save a two-dimensional mesh.
    fprintf(fout, "%d  %d  %d\n",NumTriaFace,3,trifacemarkerlist ? 1 : 0);
    for (i = 0; i < NumTriaFace; i++) {
      fprintf(fout, "%d", i + FirstNumber);
      for (j = 0; j < 3; j++) {
        fprintf(fout, "  %5d", trifacelist[i * 3 + j]);
      }
      if (trifacemarkerlist != NULL) {
        fprintf(fout, "  %d", trifacemarkerlist[i]);
      }
      fprintf(fout, "\n");
    }
  }

  fclose(fout);
}

//============================================================================//
//                                                                            //
// save_faces()    Save faces to a .face file.                                //
//                                                                            //
//============================================================================//

void tetgenio::save_faces(const char* filebasename)
{
  FILE *fout;
  char outfacefilename[FILENAMESIZE];
  int i;

  sprintf(outfacefilename, "%s.face", filebasename);
  printf("Saving faces to %s\n", outfacefilename);
  fout = fopen(outfacefilename, "w");
  fprintf(fout, "%d  %d\n", NumTriaFace, 
          trifacemarkerlist != NULL ? 1 : 0);
  for (i = 0; i < NumTriaFace; i++) {
    fprintf(fout, "%d  %5d  %5d  %5d", i + FirstNumber, trifacelist[i * 3],
            trifacelist[i * 3 + 1], trifacelist[i * 3 + 2]);
    if (trifacemarkerlist != NULL) {
      fprintf(fout, "  %d", trifacemarkerlist[i]);
    }
    fprintf(fout, "\n");
  }

  fclose(fout);
}

//============================================================================//
//                                                                            //
// save_edges()    Save egdes to a .edge file.                                //
//                                                                            //
//============================================================================//

void tetgenio::save_edges(char* filebasename)
{
  FILE *fout;
  char outedgefilename[FILENAMESIZE];
  int i;

  sprintf(outedgefilename, "%s.edge", filebasename);
  printf("Saving edges to %s\n", outedgefilename);
  fout = fopen(outedgefilename, "w");
  fprintf(fout, "%d  %d\n", NumEdge, edgemarkerlist != NULL ? 1 : 0);
  for (i = 0; i < NumEdge; i++) {
    fprintf(fout, "%d  %4d  %4d", i + FirstNumber, edgelist[i * 2],
            edgelist[i * 2 + 1]);
    if (edgemarkerlist != NULL) {
      fprintf(fout, "  %d", edgemarkerlist[i]);
    }
    fprintf(fout, "\n");
  }

  fclose(fout);
}

//============================================================================//
//                                                                            //
// save_neighbors()    Save egdes to a .neigh file.                           //
//                                                                            //
//============================================================================//

void tetgenio::save_neighbors(char* filebasename)
{
  FILE *fout;
  char outneighborfilename[FILENAMESIZE];
  int i;

  sprintf(outneighborfilename, "%s.neigh", filebasename);
  printf("Saving neighbors to %s\n", outneighborfilename);
  fout = fopen(outneighborfilename, "w");
  fprintf(fout, "%d  %d\n", NumTetr, iMeshDim + 1);
  for (i = 0; i < NumTetr; i++) {
    if (iMeshDim == 2) {
      fprintf(fout, "%d  %5d  %5d  %5d", i + FirstNumber,  neighborlist[i * 3],
              neighborlist[i * 3 + 1], neighborlist[i * 3 + 2]);
    } else {
      fprintf(fout, "%d  %5d  %5d  %5d  %5d", i + FirstNumber,
              neighborlist[i * 4], neighborlist[i * 4 + 1],
              neighborlist[i * 4 + 2], neighborlist[i * 4 + 3]);
    }
    fprintf(fout, "\n");
  }

  fclose(fout);
}

//============================================================================//
//                                                                            //
// save_poly()    Save segments or facets to a .poly file.                    //
//                                                                            //
// It only save the facets, holes and regions. No .node file is saved.        //
//                                                                            //
//============================================================================//

void tetgenio::save_poly(const char *filebasename)
{
  FILE *fout;
  facet *f;
  polygon *p;
  char outpolyfilename[FILENAMESIZE];
  int i, j, k;

  sprintf(outpolyfilename, "%s.poly", filebasename);
  printf("Saving poly to %s\n", outpolyfilename);
  fout = fopen(outpolyfilename, "w");

  // The zero indicates that the vertices are in a separate .node file.
  //   Followed by number of dimensions, number of vertex attributes,
  //   and number of boundary markers (zero or one).
  fprintf(fout, "%d  %d  %d  %d\n", 0, iMeshDim, NumPointAttribute,
          pPointMarkerList != NULL ? 1 : 0);

  // Save segments or facets.
  if (iMeshDim == 2) {
    // Number of segments, number of boundary markers (zero or one).
    fprintf(fout, "%d  %d\n", NumEdge, edgemarkerlist != NULL ? 1 : 0);
    for (i = 0; i < NumEdge; i++) {
      fprintf(fout, "%d  %4d  %4d", i + FirstNumber, edgelist[i * 2],
              edgelist[i * 2 + 1]);
      if (edgemarkerlist != NULL) {
        fprintf(fout, "  %d", edgemarkerlist[i]);
      }
      fprintf(fout, "\n");
    }
  } else {
    // Number of facets, number of boundary markers (zero or one).
    fprintf(fout, "%d  %d\n", NumFacet, pFacetMarkerList != NULL ? 1 : 0);
    for (i = 0; i < NumFacet; i++) {
      f = &(pFacetList[i]);
      fprintf(fout, "%d  %d  %d  # %d\n", f->NumPolygon,f->NumHole,
            pFacetMarkerList != NULL ? pFacetMarkerList[i] : 0, i + FirstNumber);
      // Output polygons of this facet.
      for (j = 0; j < f->NumPolygon; j++) {
        p = &(f->pPolygonList[j]);
        fprintf(fout, "%d  ", p->NumVertex);
        for (k = 0; k < p->NumVertex; k++) {
          if (((k + 1) % 10) == 0) {
            fprintf(fout, "\n  ");
          }
          fprintf(fout, "  %d", p->pVertexList[k]);
        }
        fprintf(fout, "\n");
      }
      // Output holes of this facet.
      for (j = 0; j < f->NumHole; j++) {
        fprintf(fout, "%d  %.12g  %.12g  %.12g\n", j + FirstNumber,
           f->pHoleList[j * 3], f->pHoleList[j * 3 + 1], f->pHoleList[j * 3 + 2]);
      }
    }
  }

  // Save holes.
  fprintf(fout, "%d\n", NumHole);
  for (i = 0; i < NumHole; i++) {
    // Output x, y coordinates.
    fprintf(fout, "%d  %.12g  %.12g", i + FirstNumber, pHoleList[i * iMeshDim],
            pHoleList[i * iMeshDim + 1]);
    if (iMeshDim == 3) {
      // Output z coordinate.
      fprintf(fout, "  %.12g", pHoleList[i * iMeshDim + 2]);
    }
    fprintf(fout, "\n");
  }

  // Save regions.
  fprintf(fout, "%d\n", NumRegion);
  for (i = 0; i < NumRegion; i++) {
    if (iMeshDim == 2) {
      // Output the index, x, y coordinates, attribute (region number)
      //   and maximum area constraint (maybe -1).
      fprintf(fout, "%d  %.12g  %.12g  %.12g  %.12g\n", i + FirstNumber,
              pRegionList[i * 4], pRegionList[i * 4 + 1],
              pRegionList[i * 4 + 2], pRegionList[i * 4 + 3]);
    } else {
      // Output the index, x, y, z coordinates, attribute (region number)
      //   and maximum volume constraint (maybe -1).
      fprintf(fout, "%d  %.12g  %.12g  %.12g  %.12g  %.12g\n", i + FirstNumber,
              pRegionList[i * 5], pRegionList[i * 5 + 1],
              pRegionList[i * 5 + 2], pRegionList[i * 5 + 3],
              pRegionList[i * 5 + 4]);
    }
  }

  fclose(fout);  
}

//============================================================================//
//                                                                            //
// save_faces2smesh()    Save triangular faces to a .smesh file.              //
//                                                                            //
// It only save the facets. No holes and regions. No .node file.              //
//                                                                            //
//============================================================================//

void tetgenio::save_faces2smesh(char* filebasename)
{
  FILE *fout;
  char outsmeshfilename[FILENAMESIZE];
  int i, j;

  sprintf(outsmeshfilename, "%s.smesh", filebasename);
  printf("Saving faces to %s\n", outsmeshfilename);
  fout = fopen(outsmeshfilename, "w");

  // The zero indicates that the vertices are in a separate .node file.
  //   Followed by number of dimensions, number of vertex attributes,
  //   and number of boundary markers (zero or one).
  fprintf(fout, "%d  %d  %d  %d\n", 0, iMeshDim, NumPointAttribute,
          pPointMarkerList != NULL ? 1 : 0);

  // Number of facets, number of boundary markers (zero or one).
  fprintf(fout, "%d  %d\n", NumTriaFace, 
          trifacemarkerlist != NULL ? 1 : 0);

  // Output triangular facets.
  for (i = 0; i < NumTriaFace; i++) {
    j = i * 3;
    fprintf(fout, "3  %d %d %d", trifacelist[j], trifacelist[j + 1], 
            trifacelist[j + 2]);
    if (trifacemarkerlist != NULL) {
      fprintf(fout, "  %d", trifacemarkerlist[i]);
    }
    fprintf(fout, "\n");
  }

  // No holes and regions.
  fprintf(fout, "0\n");
  fprintf(fout, "0\n");

  fclose(fout);  
}

//============================================================================//
//                                                                            //
// readline()   Read a nonempty line from a file.                             //
//                                                                            //
// A line is considered "nonempty" if it contains something more than white   //
// spaces.  If a line is considered empty, it will be dropped and the next    //
// line will be read, this process ends until reaching the end-of-file or a   //
// non-empty line.  Return NULL if it is the end-of-file, otherwise, return   //
// a pointer to the first non-whitespace character of the line.               //
//                                                                            //
//============================================================================//

char* tetgenio::readline(char *string, FILE *infile, int *linenumber)
{
  char *result;

  // Search for a non-empty line.
  do {
    result = fgets(string, INPUTLINESIZE - 1, infile);
    if (linenumber) (*linenumber)++;
    if (result == (char *) NULL) {
      return (char *) NULL;
    }
    // Skip white spaces.
    while ((*result == ' ') || (*result == '\t')) result++;
    // If it's end of line, read another line and try again.
  } while ((*result == '\0') || (*result == '\r') || (*result == '\n'));
  return result;
}

//============================================================================//
//                                                                            //
// findnextfield()   Find the next field of a string.                         //
//                                                                            //
// Jumps past the current field by searching for whitespace or a comma, then  //
// jumps past the whitespace or the comma to find the next field.             //
//                                                                            //
//============================================================================//

char* tetgenio::findnextfield(char *string)
{
  char *result;

  result = string;
  // Skip the current field.  Stop upon reaching whitespace or a comma.
  while ((*result != '\0') && (*result != ' ') &&  (*result != '\t') && 
         (*result != ',') && (*result != ';')) {
    result++;
  }
  // Now skip the whitespace or the comma, stop at anything else that looks
  //   like a character, or the end of a line. 
  while ((*result == ' ') || (*result == '\t') || (*result == ',') ||
         (*result == ';')) {
    result++;
  }
  return result;
}

//============================================================================//
//                                                                            //
// readnumberline()   Read a nonempty number line from a file.                //
//                                                                            //
// A line is considered "nonempty" if it contains something that looks like   //
// a number.  Comments (prefaced by `#') are ignored.                         //
//                                                                            //
//============================================================================//

char* tetgenio::readnumberline(char *string, FILE *infile, char *infilename)
{
  char *result;

  // Search for something that looks like a number.
  do {
    result = fgets(string, INPUTLINESIZE, infile);
    if (result == (char *) NULL) {
      return result;
    }
    // Skip anything that doesn't look like a number, a comment, 
    //   or the end of a line. 
    while ((*result != '\0') && (*result != '#')
           && (*result != '.') && (*result != '+') && (*result != '-')
           && ((*result < '0') || (*result > '9'))) {
      result++;
    }
    // If it's a comment or end of line, read another line and try again.
  } while ((*result == '#') || (*result == '\0'));
  return result;
}

//============================================================================//
//                                                                            //
// findnextnumber()   Find the next field of a number string.                 //
//                                                                            //
// Jumps past the current field by searching for whitespace or a comma, then  //
// jumps past the whitespace or the comma to find the next field that looks   //
// like a number.                                                             //
//                                                                            //
//============================================================================//

char* tetgenio::findnextnumber(char *string)
{
  char *result;

  result = string;
  // Skip the current field.  Stop upon reaching whitespace or a comma.
  while ((*result != '\0') && (*result != '#') && (*result != ' ') && 
         (*result != '\t') && (*result != ',')) {
    result++;
  }
  // Now skip the whitespace and anything else that doesn't look like a
  //   number, a comment, or the end of a line. 
  while ((*result != '\0') && (*result != '#')
         && (*result != '.') && (*result != '+') && (*result != '-')
         && ((*result < '0') || (*result > '9'))) {
    result++;
  }
  // Check for a comment (prefixed with `#').
  if (*result == '#') {
    *result = '\0';
  }
  return result;
}
